clear;
load ../Bathymetry/collect_all.mat; %Observations
load ../FLAC/spacingld.mat; %FLAC output
goffpoints=load('../Bathymetry/goff_data.csv');


%%flags 
q.diff=1;   % do computations on differenced bathymetry 
q.fmax=1/5; % max freq. to fit to
q.dt=5;
q.nw=5;
pl=find(All.S>-inf);
options=optimoptions('lsqnonlin');
options.StepTolerance=1e-20;
options.TolFun=1E-10;
options.MaxFunEvals = 1000;
options.MaxIterations = 500;

%%subset data
r=unique(All.region); 
time=All.time;
for k=1:length(r);
    pl=find(All.region==r(k));
    sr=[];
    for ct=1:length(pl);         
        z=All.z(pl(ct),:);
        z(isnan(z))=nanmean(z);
        sr(ct)=All.S(pl(ct))*1/(100*1000)*1000; %cm/year*km/cm*year/ky=km/ky        

        % Use Differences? 
        if q.diff
            y=diff(detrend(z(1:end-5)));
        else
            y=detrend(z(1:end));
        end
        y=smoothPH(y,q.dt);
        y=y(1:q.dt:end);         
        [Pw,fq]=pmtmPH(y,q.dt,q.nw,0);
        R(k).Pall(ct,:)=Pw; 
    end;
    R(k).Pavg=nanmean(R(k).Pall); 
    R(k).savg=mean(sr(ct));    
    pl2=1:length(fq); %exclude lowest freqs. as these are biased
    [dummy pl]=max(R(k).Pavg(pl2));
    R(k).F=fq(pl2(pl));
end;

%Make Figure 2 of the main text
figure(2); clf; hold on;
plot([2.3 2.3],[0 9],'k--','linewidth',1); 
plot([3.8 3.8],[0 9],'k--','linewidth',1); 
for ct=1:length(R); 
    S(ct)=R(ct).savg*100;
    Lam(ct)=R(ct).savg./R(ct).F; 
end;

%Regress across intermediate spreading ridges
pl1=[5 8:17];
[B1,BINT] = regress(Lam(pl1)',S(pl1)');
x=[0 6];
xo=3.25;
yo=xo*B1;
fill([0 6 6 0],[0 BINT(1)*x(2) BINT(2)*x(2) 0],[0.9 0.9 0.8],'edgecolor',[0.92 0.92 0.85]);
[B1 BINT];
[B1g,BINTg] = regress([Lam(pl1)'; goffpoints(6:7,2)],[S(pl1)'; goffpoints(6:7,1)]);

%Regress across fast spreading ridges
pl2=[1:4 6];
[B1,BINT] = regress(Lam(pl2)',S(pl2)');
x=[0 9];
xo=3.25;
yo=xo*B1;
fill([0 9 9 0],[0 BINT(1)*x(2) BINT(2)*x(2) 0],[0.9 0.9 0.8],'edgecolor',[0.92 0.92 0.85]);
[B1 BINT];
[B1g,BINTg] = regress([Lam(pl2)'; goffpoints(8:10,2)],[S(pl2)'; goffpoints(8:10,1)]);

%Expected from sea-level hypothesis
plot([0 9],[0 9]*41/100,'k','linewidth',1);
plot([0 6],[0 6]*100/100,'k','linewidth',1);

%FLAC results
pl=find(F.U1*100<=6);
L(1)=plot(F.U1(pl)*100,F.U1(pl).*F.pr100(pl),'s','markersize',8,'linewidth',1,'color',[0.8 0 0]);
pl=find(F.U1*100<=9);
L(2)=plot(F.U1(pl)*100,F.U1(pl).*F.pr41(pl),'bs','markersize',8,'linewidth',1);

L(3)=plot(goffpoints(:,1),goffpoints(:,2),'o','markersize',8, ...
          'linewidth',1,'markerfacecolor',[0 0.75 0.5],'markeredgecolor',[0 0.75 0.5]);

%Observed
for ct=1:length(S),
    if ct~=7,
        L(4)=plot(S(ct),Lam(ct),'o','markersize',8,'linewidth',1,'markerfacecolor',[0.5 0.5 0],'markeredgecolor',[0.5 0.5 0]);
    else,
        plot(S(ct),Lam(ct),'o','markersize',8,'linewidth',2,'markerfacecolor',[1 1 1],'markeredgecolor',[0.5 0.5 0]);
    end;
    text(S(ct),Lam(ct),['  ',num2str(ct)],'fontsize',12);
end;

set(gca,'fontsize',14);
xlabel('half-spreading rate (cm/yr)','fontsize',15); 
ylabel('characteristic wavelength (km)','fontsize',15); 
h=axis; 
h=legend(L,'FLAC 100 ky','FLAC 41 ky','Observed (Goff, 2015)','Observed (this study)','fontsize',14);

return;

%%----------------------------------------------------
%Supplemental figure: Spectogram
for ct=1:17;
    ptemp(ct,:)=R(ct).Pavg/max(R(ct).Pavg);
    wl(ct,:)=fq/R(ct).savg;
end;
wli=linspace(min(wl(:)),2,1.5*length(fq));
for ct=1:17;
    wtemp(ct,:)=interp1(wl(ct,:),ptemp(ct,:),wli);
end;

%sort according to spreading rate
Sr=[R.savg];
[dummy pls]=sort(Sr);

figure(102); clf; hold on;
subplot(3,2,[1 3]); hold on;
imagesc(wli,1:17,wtemp(pls,:));
for ct=1:17,
    [dummy pl]=max(wtemp(pls(ct),:));
    plot(wli(pl),ct,'w.','markersize',14);
end;
axis tight;
set(gca,'ydir','reverse','ytick',[1:17]);
h=ylabel('region'); 
h=xlabel('inverse wavelength (1/km)');
hc=colorbar('horiz');
hc.Label.String='normalized spectral energy';
hc.FontSize=15;
hc.Label.FontSize=17;
p=hc.Position;
p(1)=p(1)+0.2;
p(2)=p(2)-0.2;
hc.Position=p;
set(gca,'clim',[0 1]);
set(gca,'yticklabel',pls);
font(gca,15);

subplot(3,2,[2 4]); hold on;
imagesc(fq,1:17,ptemp(pls,:));
for ct=1:17,
    [dummy pl]=max(ptemp(pls(ct),:));
    plot(fq(pl),ct,'w.','markersize',14);
end;
plot([1/100 1/100],[1 17],'w--','linewidth',1.5); 
plot([1/41 1/41],[1 17],'w--','linewidth',1.5); 
plot([1/20 1/20],[1 17],'w--','linewidth',1.5); 
axis tight;
set(gca,'ydir','reverse','ytick',[1:17]);
font(gca,12);
h=ylabel('region'); 
h=xlabel('frequency (1/ky)');
%colorbar('horiz');
set(gca,'clim',[0 1]);
set(gca,'yticklabel',pls);
c=demcmap(0);
c(1:32,:)=sqrt(sqrt(c(1:32,:)));
colormap(c);
font(gca,15);

